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NUMERICAL  STUDIES  OP  FRONTAL  MOTION 
IN  THE  ATMOSPHERE   -   I 

ABSTRACT 

The  motion  of  frontal  disturbances  in  the 
atmosphere  is  studied  by  the  numerical  solution  of 
differential  equations  based  upon  a  two-layer  model 
of  an  incompressible  fluid  on  a  rotating  earth.   The 
density  of  each  layer  is  assumed  to  be  constant.   The 
upper  and  lower  fluids  correspond  respectively  to  warm 
and  cold  air.   Tn  this  first  attempt,  only  the  motion 
of  the  lower  cold  air  layer  is  studied  by  assomiing,  in 
effect,  that  the  dynamics  of  the  perturbations  in  the 
upper  warm  air  layer  can  be  neglected. 

The  numerical  study  of  this  simple  mechanical  model 
shows  that  even  though  thermodynamic  processes  have  been 
ignored,  the  occlusion  process,  characteristic  for  warm 
and  cold  fronts,  develops  from  an  initially  sinusoidal 
frontal  pattern.   Two  cases  of  different  initial 
conditions  are  examined.   Case  A:   Only  the  east-west 
component  of  wind  velocity  is  initially  geostrophlc. 
Case  B:   Both  east-west  and  north-south  components  are 
initially  geostrophlc.   In  both  cases,  computations 
indicate  that  the  cold  front  propagates  faster  than  the 
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warm  front  and  that  a  relatively  strong  mass  convergence 
zone  appears  behind  the  cold  front  only.   This  fact  suggests 
the  occurrence  of  severe  storms  associated  with  cold  fronts, 
but  not  with  warm  fronts  in  the  atmosphere.   The  numerical 
method  developed  here  to  calculate  the  movement  of  the 
front  is  based  on  following  the  motion  of  the  material 
"particles"  at  the  front.   This  method  has  applications 
to  the  numerical  solution  of  a  certain  class  of  hydro- 
dynamic  flow  problems  in  which  the  entire  boundary  of  the 
domain  of  integration  is  not  given  a  priori,  but  must  be 
determined  (so-called  free-boundary  problems). 
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NUMERICAL  STUDIES  OP  FRONTAL  MOTION 

IN  THE  ATMOSPHERE -I 

1.   Introduction. 

It  is  well  known  in  meteorology  that  the  weather  in 
middle  latitudes  is  determined  primarily  by  a  series  of 
events  associated  with  the  development  and  propagation  of 
wave-like  disturbances  (frontal  cyclones)  on  a  disconti- 
nuity surface  ("polar  front")  between  warm  (tropical)  and 
cold  (polar)  air  masses  in  the  atmosphere  (Bjerknes,  J. 
and  H.  Solberg,  1922) . 

The  theoretical  study  of  frontal  cyclones  was  first 
undertaken  by  Solberg  (I928)  using  the  method  of  small 
perturbations  (see,  also  V.  Bjerknes  et  al.,  1953)'   The 
main  results  of  Solberg 's  theory  which  are  relevant  to  the 
cyclone  problem,  were  recapitulated  by  J.  Bjerknes  and 
Gods-ke  (I936).   Mathematical  formulation  of  the  theory 
was  improved  later  by  Kotschin  (1932)  and  recently  by 
Eliasen  (I960)  by  taking  into  account  the  correct  boundary 
conditions  for  the  problem.  The  aim.  of  these  linearized 
studies  was  to  find  solutions  whose  patterns  resemble  those 

The  research  reported  in  this  paper  was  performed  under 
Contract  Nonr-285(55)  with  the  Office  of  Naval  Research, 
U.  S.  Navy  and  under  Contract  AT(30-l)-l480  with  the 
U,  S.  Atomic  Energy  Commission.   Reproduction  in  whole 
or  in  part  is  permitted  for  any  purpose  of  the 
United  States  Government. 
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of  observed  nascent  cyclones c   Since  the  initial  phase 
of  a  frontal  cyclone  could  be  identified  with  a  wave-like 
disturbance  of  small  amplitude  superimposed  upon  a  more  or 
less  straight  quasi-stationary  front,  it  was  natural  to 
begin  by  considering  motions  which  depart  so  little  from 
such  a  steady  motion  that  linearizations  of  the  relevant 
equations  can  be  performed.   However,  this  linearization 
limits  the  study  of  the  complete  evolution  of  a  develop- 
ing cyclone  since  it  holds  only  for  small  departures  from 
the  steady  motion. 

-  The  im.portance  of  nonlinear  effects  in  the  dynamical 
equations  which  describe  the  movement  of  a  cold  front  were 
pointed  out  by  Freeman  (1952),  Abdullah  (19^9)  and  Tapper 
(1952)  who  applied  the  method  of  characteristics  to  solve 
relevant  nonlinear  equations  with  one  space  variable  and 
time  (one-dimensional  problems) »   Present  computing 
machines  and  difference  methods  for  solving  multi-dimensional 
hydrodynamic  flow  problems  enable  us  to  integrate  more 
complicated  nonlinear  equations  numerically.   With  this 
stiuation  in  mind.  Stoker  (1553; 7)  developed  a  two-layer 
model  for  the  study  of  the  motion  of  a  frontal  surface 
(see.  Section  2).   Whitham  (1953)  niade  a  qualitative  study 
of  this  model  which  indicated  strongly  that  the  evolution 
of  frontal  cyclones  might  well  follow  the  pattern  that  leads 
to  the  occlusion  process,, 

The  object  of  this  paper  is  to  present  a  finite-difference 
method  for  dealing  with  this  model  and  to  discuss  the  results 


of  numerical  integration  under  prescribed  Initial  and 
boundary  conditions.   In  doing  so,  the  model  was  further 
simplified  by  assuming,  in  effect,  that  the  dynamics  of 
the  perturbations  in  the  warm  air  layer  could  be  neglected. 
With  this  assumption  the  differential  equations  reduce  to 
a  system  of  three  equations  for  three  dependent  variables 

* 

(two  horizontal  velocity  components  and  the  height  of 
the  cold  air  layer)  with  three  independent  variables  (two 
space  coordinates  and  tim.e)  as  discussed  in  Section  2.   A 
special  numerical  difficulty  arises  from  the  fact  that  the 
front  is  a  free  boundary  along  which  the  differential 
equations  are  in  a  certain  sense  singular.   The  first 
numerical  calculations  based  on  this  formulation  were  made 
by  Fife  (1959)  who  considered  various  initial  conditions  and 
boundary  conditions  for  the  problem,  except  for  the  conditions 
at  the  front  along  which  nothing  was  prescribed;  it  was 
fixed  simply  as  the  line  along  which  the  depth  of  the  cold 
air  layer  is  zero.   The  limitations  of  the  Univac  I  computer 
did  not  permit  Fife  to  reduce  the  interval  size  sufficiently 
to  complete  the  numerical  study.   This  problem  was  overcome 
by  Lewis  (I961)  who  used  the  IBM  709O  computer,   Lewis 
improved  the  evaluation  of  derivatives  at  the  frontal 
boundary  with  a  method  which  is  accurate  to  first  order  in 
the  space  increment.  Results  of  his  calculations,  however, 
indicated  the  need  of  a  more  accurate  procedure  at  the 
boundary.   The  scheme  presented  in  this  paper  was  therefore 
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devised  so  that  it  is  correct  to  second  order  in  the  space 
increment  at  the  boundary  as  well  as  in  the  interior  of  the 
domain.   The  position  of  the  frontal  boundary  is  calculated 
following  the  movement  of  the  material  "particles"  cons_^titut- 
ing  the  boundary.   The  numerical  calculations  were  carried 
out  over  an  interval  of  time  large  enough  so  that  the 
asymmetric  shape,  expected  from  actual  observations  of 
frontal  motions,  developed  from  an  initially  sinusoidal 
pattern. 

2.   Differential  equations  for  frontal  motion. 

Figure  la  shows  the  initial  state  of  the  dynamical 
system  to  be  studied  here.   It  shows  a  cold  wedge  of  air 
at  the  ground  with  a  warm,  layer  over  it .   The  fluid  in 
each  layer  is  assumed  to  be  a  perfect  incompressible  fluid 
with  constant  density  subject  to  gravity.   Initially  the 
velocity  in  both  layers  is  constant  and  in  the  direction 
of  the  X-axis  (to  be  thought  of  as  the  eastward  direction), 
but  it  will  in  general  have  different  values  in  the  two 
layers  so  that  the  interface  between  the  layers  is  not  only 
a  discontinuity  surface  for  the  density,  but  is  also  a 
surface  across  which  a  discontinuity  in  the  tangential 
velocity  components  in  the  two  layers  may  occur.   The 
velocity  in  the  warm  air  layer  would  be  that  of  the  pre- 
vailing westerlies;  in  the  cold  air  the  m.otion  might  be  in 
the  opposite  direction.   The  reason  that  the  interface  is 
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Initially  a  plane  Inclined  to  the  horizontal  is  that 
the  coordinate  system  is  assumed  to  be  rotating  about 
the  vertical  axis  with  a  constant  angular  velocity  O. 
In  effect,  the  rotation  of  the  earth  is  taken  into  account 
by  assuming  that  O  =  cd  sin  ^     =     -p  f   (f,  the  Coriolis 
parameter),  with  oj  the  angular  velocity  of  the  earth  and 
(})  the  latitude  of  the  origin  of  the  coordinate  system; 
this  means  that  the  sphericity  of  the  earth  is  not  taken 
into  account. 

In  addition  to  the  assumptions  already  made.  I.e. 
that  the  cold  and  warm  layers  consist  of  perfect,  incom- 
pressible, gravitating  fluids,  a  basic  assumption  commonly 
made  in  m.eteorology  in  discussing  large  scale  motions  is 
also  made  here,  i.e.  that  the  hydrostatic  pressure  law 
holds.   This  is  equivalent  to  assuming  that  the  vertical 
component  of  the  acceleration  of  the  fluid  particles  is 
negligible.   With  reference  to  Figure  lb,  which  shows  a 
north-south  vertical  section  of  the  fluid,  the  pressure 
in  the  two  layers  is  given  by  the  formulas 

p'(x,y,z,t)   =   p'g(h»  -  z) 

(2.1) 

p(x,y,z,t)    =   p'g(h'-h)  +  pg(h-z)  . 

As  in  all  that  follows,  all  analogous  quantities,  such  as 
the  pressure  p,  density  p,  height  h  of  the  upper  surface 
of  a  layer,  velocity  components  u,  v,  are  distinguished  in 
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the  two  layers  by  using  primes  for  those  in  the  warm  layer. 

Since  both  h  and  h'  are  functions  of  x,  y,  and  t  only, 

while  p,  p'  and  g  are  constants,  it  follows  that  the 

horizontal  pressure  gradient  (p^,p^J  depends  only  upon 

X  y 

X,  y,  and  t  and  not  on  z.   Since  the  Coriolis  force  may 
also  be  regarded  without  serious  error  as  being  independent 
of  z  it  follows  that  the  horizontal  acceleration  of  all 
particles  is  independent  of  z,  and  hence  that  the  horizontal 
components  will  remain  independent  of  z  if  that  were  to  be 
true  initially  —  which  we  assume  here.   Hence  u  -   u(x,y,t), 
V  =  v(x,y,t)  and  the  Euler  equations  of  motion  in  the  two 
layers  are  as  follows,  when  the  pressure  gradients  are 
computed  from  (2.1): 


(2.2) 


(2.3) 


u,  +  uu  +  vu   =   -  g[-^  h   +  (1  -  -£-h  1  +  fv 
t     X     y       ^'■p   X    ^     p   X-' 


V,  +  uv  +  vv   =   -  g[£-  h   +(l--£-h]-fu 
t     X     y       ^'•p   y   '    P   y 


u'  +  u'u*  +  v'u'   =   -  gh  +  fv' 
t      X      y       ^  X 


I      II      It  t       r 

v^_+uv^  +  vv   =   -gh   -fu 
t      X      y       ^  y 


Since  the  fluids  in  both  layers  are  assumed  to  be 
incompressible,  the  equations  of  continuity  to  be  added 
to  the  dynamical  equations  are 
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(2.^) 


(uh)^  +  (vh)   +  h   =  0   , 


[u'(h'-h)]^  +  [v'(h'-h)]^^  +  (h'-h),   =  0 


'X 


y 


t 


These  equations  simply  state  that  the  flux  of  fluid  Into  a 
vertical  column  of  the  fluid  Is  just  balanced  by  the  change 
of  height  of  the  column. 

It  is  Important  to  point  out  that  this  system  of 
partial  differential  equations  has  an  exact  solution  that 
corresponds  to  the  initial  state  of  the  stationary  front 
with  which  the  discussion  has  started.   It  is  the  solution 


_  I 


u 


u    (=  constant) 


(2.5a) 


and 


v'   =   0 


_  I 

h'   =   h   ,   where 


ah'/^x  =  0 


V  Sh'/Sy 


-  f  u' 


u  =   u   (=  constant) 


V  =  0 


(2.5b)  \   h  =   h   ,   where 


Bh/Sx  =  0 


I  _  1 


V  afi/3y  =  g(i4vp)  'f  =  -  "' 


Note  that  even  though  p '/p  is  close  to  1  in  value  since  the 
cold  and  warm  layers  differ  in  teiiiperature  by  only  a  few 
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degrees.  It  nevertheless  turns  out  that  h   is  small  since 
the  Coriolis  parameter  f  is  very  sm.all  compared  with  g, 
and  slopes  a  of  the  discontinuity  surface  are  of  the  order 
of  a  fraction  of  a  degree  when  parameter  values  correspond- 
ing to  those  of  the  atmosphere  are  chosen  (see.  Section  8). 

The  last  fact,  in  combination  with  some  of  the  basic 
assumptions  made  here,  makes  it  seem  reasonable  to  simplify 
the  mathematical  problem  still  more  without  necessarily 
impairing  the  approximation  to  the  m.otion  of  the  system. 
The  system  of  differential  equations  (2.2),  (2.3),  (2.4), 
together  with  appropriate  Initial  and  boundary  conditions, 
is  now  under  investigation  by  means  of  numerical  methods. 
As  one  sees,  it  is  a  system  of  nonlinear  partial  differen- 
tial equations  for  six  functions  u,  u',  v,  v',  h,  h'  depend- 
ing on  three  variables  x,  y,  and  t.   This  paper,  however,  is 
concerned  with  an  approximate  treatment  that  cuts  the  number 
of  dependent  variables  in  half  by  simply  neglecting  the 
dynamics  of  the  perturbations  in  the  warm  air  layer  relative 
to  the  initial  stationary  state.   This  assumption,  at  first 
sight  rather  drastic,  has  a  reasonable  physical  basis,  as 
follows.   Imagine  the  stationary  front  to  have  developed  a 
bulge  in  the  northward  direction,  say,  as  in  Figure  2a.   The 
warm  layer  can  adjust  itself  to  the  new  condition  simply 
through  a  slight  change  in  the  vertical  velocity  component 
without  a  change  in  the  horizontal  velocity  components  u' 
and  v',  as  indicated  in  Figure  2b,  which  is  a  vertical 


section  taken  along  the  line  AB.   However,  in  the  cold  air 
layer  quite  large  changes  in  u  and  v  will  certainly  be 
necessary  if  the  bulge  develops  at  all  in  the  way  it  does 
in  the  occlusion  process.   Thus  in  this  approximate  theory 
it  is  assumed  that  u',  v',  and  h'  retain  for  all  time  the 
values  they  had  in  the  initial  steady  state  (2.5a).   The 
differential  equations  for  u,  v,  and  h  are  therefore  the 
following: 

(2.6)  u^  +  uu^  +  ^y  +  S^l  "  f~^\     "   ^^ 

(2.7)  v^  +  uv^  +  vvy  +  g(l  -  ^)hy  =   f(£_L  u'-u) 

(2.8)  h,  +  (uh)^  +  (vh)^^  =      0   . 

i;        A       y 


3.   Boundary  and  initial  conditions. 

Our  task  is  to  solve  equations  (2.6)  and  (2.8)  with 
prescribed  boundary  and  initial  conditions  in  a  region 
given  by  (x,y)  e  D  oR,   t  >  0,  where  R  is  a  fixed 
rectangular  domain  (see,  Figure  3)  with  sides  parallel 
to  the  axes.   We  assume  the  following  conditions: 

(a)  The  cold  air  lies  over  a  region  Da   R  bounded 
by  the  curve  C  on  the  south  as  illustrated  in 
Figure  5;  the  warm  air  is  above  the  cold  air  and 
over  the  domain  D',  the  complement  of  D  in  R . 

(b)  The  curve  C  is  the  surface  front  which  we,  for 


brevity,  call  "front." 

(c)   The  curve  C  has  a  continuous  tangent. 

The  numerical  solutions  of  (2.6)  to  (2.8)  are  confined 
to  D  which  we  call  the  cold  air  mass  domain.   We  refer  to  D' 
as  the  warm  air  mass  domain.   Since  C  is  defined  as  the 
line  along  which  h  =  0  it  is  seen  from  (2.8)  that  the  line 
C  must  move  with  the  particle  velocity  [u  ,v  ]  on  C .   Let 
[x  ,y  ]  be  the  coordinates  of  a  point  on  C,  then  the 
boundary  conditions  along  C  are  that 

h  =   0 
(3.1)   (  ^    (x^(t))   =  u^(x^,ye^t) 

It  (^c^^))  -     ^c^^c'^c^^) 

where  d/dt  signifies  the  material  derivative. 

The  boundary  conditions  on  the  other  three  boundaries 
for  u,  V,  and  h  could  be  prescribed  in  a  variety  of  ways, 
depending  on  the  physical  conditions  assumed  to  hold  there. 
For  the  sake  of  simplicity  in  making  numerical  calculations, 
we  decided  to  assume  throughout  this  paper  that 

(1)  The  y-component  of  the  velocity  vanishes  for 
all  time  along  the  northern  boundary. 

(2)  u,  V  and  h  are  periodic  in  the  space  variable 

X  with  a  period  equal  to  the  distance  between 

^According  to  the  Glossary  of  Meteorology  (American  Meteoro- 
logical Society,  1959),  the  term  "front"  is  applied  to  the 
interface  or  transition  zone  between  two  air  masses  of 
different  density.   The  line  of  intersection  of  the  disconti- 
nuity surface  with  the  earth's  surface  is  called  a  "surface 
front".  However,  in  this  paper,  we  follow  a  customary  usage 
and  "surface  front"  is  just  "front." 
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the  east  and  west  boundaries. 

Since  the  atmospheric  motion  on  the  sphere  is  periodic  in 
the  longitude,  condition  (2)  isn't  too  far  removed  from 
meteorological  significance.   Furthermore,  if  the 
"occlusion"  takes  place  in  a  relatively  small  region 
centered  in  the  period,  we  may  by  varying  the  size  of  the 
period  determine  the  influence  of  the  periodicity  assump- 
tion.  The  effect  of  varying  the  boundary  conditions  upon 
the  motion  of  the  cold  air  remains  to  be  investigated. 

In  order  to  solve  an  initial  value  problem  for 
equations  (2.6)  to  (2.8),  the  values  of  u,  v  and  h  in 

the  cold  air  domain  D  and  u  and  v  on  the  front  C  must 

c      c 

be  given  at  the  initial  time  t  =  0.   In  this  paper,  the 
following  two  sets  of  initial  conditions  are  considered. 
These  initial  conditions  may  be  regarded  as  perturbations 
on  the  steady  state  described  by  (3.2). 

(1)  Case  A:   Only  the  east-west  component  of  the  cold 
air  velocity  is  initially  geostrophic.   The  initial 
conditions  are  that 

u  =  u   (=  constant) 

V  =   0 


f 


^^•2)     h  =  gd^yp)  (-^u  -u)  (y-y^),  f or  Y  >  y  > 


y. 


in  which 
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-£—  u   >  u   ,        (u'  =  constant) 

r 

y^  =   C^  sm  (—X-  -^)    +  C^      , 


where  X  and  Y  denote  the  lengths  of  the  domain  R  In  the 
X  and  y-dlrections .   Note  that  the  slope  of  h  is  constant 
in  the  y-direction,  but  varies  sinusoidally  in  the 
x-direction.   The  equation  for  y  determines  the  y 
coordinate  of  the  front  as  a  function  of  x;  C^  and  Cp 
are  parameters. 

(2)   Case  B:   Both  the  x  and  y  components  of  velocity 
are  initially  geostrophic. 

We  assume  the  following  form  for  the  height  of  the 
cold  air  layer  at  t  =  0 . 

y-y 

h  =  {-.T-rr)^  >  f or  Y  >  y  >  y,  ,  with 

(3.3a)    H  =   g^fp.yp^  (^  n'-u)  (Y-b)  ,  £l  u'  >  u,  and 

y   as  given  in  (3-2) . 

Here  H  denotes  the  height  of  cold  air  at  the  northern 
boundary.   The  height  h  is  linear  in  y,  but  the  slope 
^h/Sy  is  a  function  of  the  x-coordinate .   The  x  and  y 
components  of  the  cold  air  velocity  at  t  =  0  are  computed 
from 

The  use  of  an  initially  geostrophic  wind  field  was 
suggested  by  Norm.an  Phillips. 
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(3.3b) 


^  _  g(l-p'/p)  4h  ^21   u'  , 


V  _   g(l-p'/p)  ^h 
^  "       f     ^ 


This  Initial  wind  field  is  geostrophic  in  the  sense  that 
substitution  of  (3.3b)  in  (2.6)  and  (2.7)  gives  du/dt  =  0 
and  dv/dt  =  0  where  d/dt  signifies  the  particle  derivative. 

4.   Energy  equation  for  the  cold  air. 

Upon  multiplying  (2.3)  by  p  and  integrating  the  result- 
ing equation  with  respect  to  x  and  y  over  the  domain  D  while 
invoking  the  boundary  conditions  as  discussed  in  Section  J), 
one  obtains  the  conservation  equation  for  the  total  cold 
air  mass,  M 

(4.1)  1^  =   0    with   M  =   p  Th  dS 

D 
where  dS  =  dx  dy. 

Likewise,  upon  multiplying  (2.6),  (2.7)  and  (2.8) 
by  phu,  phv  and  gp(l-  -^)h,  respectively,  and  then  integrat- 
ing the  resulting  equations  over  D  and  adding  them  together, 
one  finds  that 

(4.2)  -^    (K+P)   =  W 
where 


K  =  I  y  (u^+  v^)h  dS 


D 
is  the  total  kinetic  energy  of  the  cold  air; 

-I 
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I   g(l-  £^)  J\2   ^3 


is    the   total  potential   energy  of   the   cold  air  and 

W     =      fp'u'/vhdS 

denotes  the  rate  of  work  done  upon  the  cold  layer  by  the 
action  of  the  Corlolis  force  in  the  warm  layer.   Thus, 
(4.2)  states  that  the  sum  of  the  total  kinetic  and 
potential  energy  of  the  cold  air  changes  at  a  rate 
proportional  to  the  correlation  between  the  v  and  h 
patterns  of  the  cold  air.   We  will  show  later,  however, 
that  the  magnitude  of  W  is  very  small  for  the  initial 
conditions  applied  here,  and  that  the  sum  of  P  and  K 
remained  approximately  constant. 

5.   Moving  coordinates  and  dimensionless  variables . 

For  the  purpose  of  programming  the  present  model  for 
high-speed  computers,  it  is  convenient  to  transform  the 
set  of  equations  (2.6)  to  (2.8)  by  using  the  following 
new  dimensionless  variables, 

T  =   t/At 

^   =   (x-  ut)/As  ,     r\      =      y/As 

(5.1)  u  =   (u-u)/(As/At)  ,  V  =   v/(As/At) 


./v  n  '   At  2 

h  =  g(l-  ^^)  {-J-)      h  ,   and  the  parameter 

,      At 

^  =   Ai 


Ik   - 


where  u  is  a  constant  having  the  dimensions  of  velocity. 
At  and  As  denote  units  for  time  and  length  with  ratio  A, 
a  suitably  determined  constant.   Also  we  define  the 
following  new  parameters: 

(5.2)  ^  -  ^  H  (^u'-u)  ,    F  ^   fAt  . 

With  the  aid  of  (5.I)  and  (5.2),  the  new  transformed 
equations  of  (2.6)- (2.8)  become 

(5.3)  u     +  uUj..    +  vu     +  h^      =      Pv 

(5.4)  V     +  uv.    +vv     +h        =      -Fu+G 

(5.5)  h      +h(u^+v)+  uh.    +  vh        =      0       . 

Equation  (5.5)  may  be  written  in  the  following  form  after 
being  multiplied  by  l/y  h  , 

(5. 50         i^  +|(u^  +  v^)  +  ^(j)^  +  v(|)^   =  0   , 

where  ^  =   2')/ h  .   The  (^,ri)  system  progresses  with  the 
velocity  u  along  the  positive  x-axis  which  coincides  with 
the  ^-axis,  while  the  rj-axis  remains  parallel  to  the  y-axis 

The  boundary  conditions  (3-1),  at  the  front,  when 
written  in  the  "moving  coordinates"  become 

dX 

(5.6)  d^  =  ^c  '     ^^"""S  ^   (h  =  0  )  , 

where 
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(5.7) 


X 


V. 


A 

u 


A 

V, 


We  may  combine  (5.3)  and  (5-4)  in  the  form; 
dV.    ^         A 


(5.8) 


where 


dT    "^  c  2    ' 


7  = 


0 


-F 


+F 


0 


K^  = 


0 


along  C  (h  =  0), 


A 

Vh 


A. 

ah 


^ 


Sh 


^ 


In  the  rest  of  this  paper,  we  deal  only  with  the 
"moving  system"  unless  otherwise  mentlonedj  and  we  shall 
omit  writing  circumflex  symbols  for  dimensionless  variables 
whenever  references  are  made  to  (5O/  -  (5->8}.   The 
prediction  of  the  flow  in  the  domain  D  is  based  upon  the 
system  of  equations  (5.3^;.  (5°^)  and  (5o')  which  may  be 
written  in  the  following  compact  form: 


(5-9) 


P   +  AP.  +  BF 

T        ^        T) 


where  P  and  K  are  vectors. 


(5.10) 


u 

V 


QP  +  K 


r 


K  =   a 
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while  A  and  B  are  symmetric  matrices  and  Q,  is  a  constant 
matrix  defined  by 


(5.11)    A  = 


u 

0 

1? 

0 

u 

0 

L# 

0 

u 

B  = 


V 

0 

0 

0 

V 

¥ 

0 

i* 

V 

0    P    0 

•F    0    0 

0    0    0 


The  system  (5.9)  is  hyperbolic  with  the  characteristic 
cone  consisting  of  the  instantaneous  particle  path 

(C-Cq)  =  '^0^'^"'^0^'   (^-^o)  "   ^0^'"^~'^0^  ^^'^   ^^®  circular 
cone  of  radius  Y^(T-Tp^)   centered  about  the  instantaneous 

particle  path  (Courant  and  Hilbert,  Vol.  II,  I962) .   Thus 

any  surface  tangent  at  each  point  {i,r],x)    to  either  a 

particle  path  or  to  the  circular  cone  described  above  will 

be  characteristic.   Of  course,  at  the  front  where  h  =  0, 

the  two  types  of  characteristic  surface  coalesce  since 

the  circular  cone  degenerates  into  the  particle  path. 

6.   Difference  equations  for  regular  net  points. 

Consider  the  domain  R  with  sides  of  length  L-j^  and  L^ 
in  the  {^  ,r])    system..   We  choose  a  rectangular  grid  R.  in 
such  a  way  that  the  coordinates  of  the  grid  points  are 
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defined  by  i-    =   iL-,/(IAs)^   t]  .  =    JLp/(JAs)   (1=0,1,..., I; 
j  =  0,1,..., J)  with  the  boundary  lines  1  -   0,1   and  J  =  0,J. 
For  actual  calculations  we  chose  the  grid  Intervals  L-,/1 
and  Lp/J  to  be  equal  to  As  so  that  ^.  and  r\  .  are  Integers. 

We  define  D.  as  a  connected  domain  of  net  points  in 
the  Interior  of  D  (I.e.  we  exclude  the  points  on  C).  By  a 
connected  domain  of  net  points,  we  mean  a  set  of  points  such 
that  each  point  In  the  set  has  at  least  one  neighboring  point 
in  the  |  or  r)  directions  belonging  to  the  set.   By  a  regular 
net  point  in  D.  we  mean  a  point  whose  eight  nearest  neighbors 
are  in  D..   Points  of  D.  that  are  not  regular  are  called 
Irregular .   Since  we  assume  periodic  boundary  conditions  at 
the  eastern  and  western  boundaries  of  D. ,  we  shall  regard 
the  grid  points  on  these  boundaries  as  regular  net  points 
in  D. .   Special  difference  equations  at  the  northern 
boundary  and  at  the  front  are  discussed  separately  in 
Section  7» 

For  the  regular  net  points  of  D.,  we  use  a  difference 
scheme  of  second  order  accuracy  in  As.   (This  scheme  would 
be  equivalent  to  the  Lax-Wendroff  (I960,  I962)    schem.e  if 
the  coefficients  were  constant.)   We  choose  the  time 
increment  to  be  At,  so  that  the  dlmenslonless  tim.e  increment 
At  is  unity.   With  A| ,  A'q  denoting  dlmenslonless  grid 
distances  (they  are  both  unity  by  cur  choice  of  the  square 
grid  size),  the  difference  equation  for  (5.9)  may  be 
derived  by  first  setting 
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riP  1 

-  <k>   II  {i,r],r   +-i  At)  At 

-  <B>^  (e.Tl.T  +^  At)  At 
+  (Q  <P>  +  K)  At 


where 

<  A>  = 


-|  fA(^,Ti,T  +  At)  +  A(^,r|,T)] 


and  similarly  for<B>  and  <P>  .   By  expanding  hP/bi 
and  SP/Sr|  at  {^  ,r[,x  +  ^   At)  into  a  Taylor  series  centered 
at  {i,r\,x)    up  to  terms  which  ensure  that  (6.1)  is  correct 
to  third  order   in  At,  we  obtain 

P(^,il,T  +  At)   =   P(4,ti,t) 


)  P    SA  SP      hi 


.  <A>  [||  -  ^  (A  n  +  1^  H  +  B  IS.^11  H  -Q  l|)]AT 


SP    At  ,„  a%  ,   hk   bP    ,  „  S^P  ,  ^B  SP   ^  ap. 


(6.2) 

-<B>[^-i|i  (A^  +  ^^+B|-|-f^  --Q^)]At 

+   (Q  <P>  +  K)At  . 

The  quantities  in  square  brackets  are  evaluated  at  (^^tj^t). 
We  replace  these  derivatives  in  (6.2)  at  the  regular  net 
points  by  using  the  following  central  difference  approxima- 
tions for  the  derivatives  

_ 
We  say  that  such  a  scheme  has  second  order  accuracy,  since 

^-  -^   X,   A      ^      SP  /■   ,  At,    P(t4-At)-P(t)  ,^,  .  2, 
if  we  divide  by  At,  and  use  -^    (t+-^)  =  — ^ -^ ^ — i-+0(AT  ) 

then  the  difference  quotient  form  of  the  equation  (6.1) 

has  second  order  accuracy. 
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^  =   A?  =   2Ar  ^-^C   "  '^^  '  '  (C  =  ^  or  Ti) 

2        2 

(6.3)  -^  ^  ^  =  — i^  (tI"^  -  2  +  T-^),(C  =  e  or  Ti) 

2        2 

^C^   ~   AlA^   "   4A^Ati  ^^e   "  ^      n   "  S  ^ 

where  T> ,  T  denote  translations  of  the  Independent  variable 
by  the  amount  sA^  and  sAri  in  the  ^  and  v^    directions 
respectively.   If  the  right  hand  side  of  (5-9)  is  Ignored, 
while  a/a  and  b/A  are  assumed  to  be  constant,  then  the 

above  difference  scheme  associated  with  (6.2)  Is  convergent 

2      2 
as  At  -and  As  ->  0  If  the  eigenvalues  of  A  and  B  are  less 

than  1/8.    We  therefore  e.xpect,  since  A  and  B  are 

proportional  to  \,    that  If  ?\  =  At/As  Is  kept  fixed  and 

small  enough,  the  dlmenslonless  velocity  and  height  will 

be  small  enough  in  0  _<  t  _<  T  to  ensure  that  the  convergence 

criterion  is  satisfied,  and  hence  that  the  difference  method 

will  converge  for  the  nonlinear  case  (6.2). 

The  fact  that  the  elem.ents  of  A  and  B  are  linear 

functions  of  the  components  of  P  permits  one  to  approximate 

(6.2)  in  the  following  explicit  form: 

(6.4)  P(4,n^^  +  At)   =   S  P  (^,ti,t'  +  K*  At 

where  S  is  a  difference  operator  of  second  order  with  respect 

to  4  and  T\,    and  K   is  a  vector.   The  form.s  of  S  and  K  are 
_ 

This  sufficient  condition  for  convergence  was  established 

by  Lax  and  Wendroff  (I962). 
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given  in  the  Appendix. 

7 •   Finite-difference  approximations  for  the  frontal 
and  northern  boundary  conditions . 

The  computations  for  solving  (6.2)  near  the  free 
boundary  (front)  and  at  the  northern  boundary  require 
special  handling.   In  conformity  with  the  use  of  a 
scheme  accurate  to  second  order  for  computations  at  the 
regular  net  points,  we  should  secure  at  least  second  order 
accuracy  in  evaluating  first  order  space  derivatives  at 
irregular  net  points.   Promi  (6.2)  we  see  that  the 
second  order  derivatives  have  the  factor  (At)   and  hence 
need  only  be  correctly  approximated  to  first  order  in  As. 
We  shall  first  discuss  a  m.ethod  of  predicting  the  front 
movement.   This  method  is,  in  principle,  similar  to  the 
one  used  by  Richtmyer  (I961)  for  shock  wave  calculations. 

We  represent  the  front  C  with  the  help  of  a  string 
C.  of  "front  points"  spaced  a  little  closer  than  the  net 
points,  as  shown  in  Figure  4„   We  number  the  front  points 
in  C.  consecutively  from  left  to  right  and  refer  to  a 
point  by  its  number.   Between  two  consecutive  front  points, 
the  front  line  C  is  represented  by  a  cubic  (that  is,  the 
front  points  between  the  points  I    and  (^+1)  lie  on  a 
cubic  through  the  points  with  indices  (i-1),  I,    (^+1), 
and  (i+2)  in  the  form  ii  =  f(0  o^  4  =  s(il)  depending  on 
the  slope).   At  front  point  I    the  coordinates  are  ^^(i^)j 
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ti.(t)  or  according  to  (5.7)  simply  X„j  and  the  particle 
velocity  is  (u^Ct),  v. (t))  or  simply  V. (t). 

In  order  to  solve  the  simultaneous  first-order 
ordinary  differential  equations  (5.6)  and  (5«8)  A'ith 
second-order  accuracy,  the  modified  m.ethod  of  Euler 
(trapezoidal  formula)  is  used.   The  difference  analogs 
of  (5.6)  and  (5.8)  are  then 


(7-1)  X_^(t+At)      =     Xj^ix)    +AT<V_g> 


(7.2)  V^(t   +  At)      =      V^(t)    +^7  <V^>   -  <Vh^>   +  K2"^At 

■where   the   symbol   <  >   applied   to  any   function  W  of  t 
signifies 

(7.3)  <W>    =     \    (W(t   +  At)    +  W(t))     . 

With  the  aid  of  notation  (7-3) p  (7-2)  may  be  written  in 
the  following  form: 


(7.4) 


V^(t    +  At) 


aV^(T)    -   P<Vhp>    +  PK^ 


where 


a     = 


P      = 


1   +(|  At)^ 


1   +(|  At)^ 


l-(- 


FAt  ^  2 
2    '' 


P'At 


.     ,FAt,2 
-FAt        j.-(-^) 


At 


F(At) 


2- 


F(At) 


At 
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Note  that  a  is  an  orthogonal  matrix  and  hence  both  eigen- 
values of  a  have  absolute  value  1  for  any  choice  of  the 
value  of  FAt .   The  two  formulae  (7.1)  and  (7.4)  are 
implicit  equations  for  X. (t+At)  and  V„(t+At).   The 
following  iteration  method  is  applied  to  solve  these 
equations . 

Step  1.   Set  Vh„(T)  for  Vh„ (t+At) . 

Step  2.   Predict  V.(t+At)  using  [J  A) . 

Step  3.   Predict  X. (t+At)  using  (7.1). 

(7.5)   Step  4.   Calculate  Vh^(T+AT)  by  a  scheme  to  be 

described  later. 
Step  5.   Repeat  steps  2,    3  and  4^  if  necessary,  ■•■•';■: 
until  no  further  changes  occur.   (¥e  found 
that  repeating  the  process  just  once  is 
sufficient  since  very  little  change  occurred 
in  the  second  repetition.) 

Now,  we  discuss  a  method  of  approximating  derivatives  of  a 
function  f  =  f{^,r\)    so  that  (6.2)  and  (7-4)  have  second-order 
accuracy  at  irregular  net  points  (i.e.  points  which  lack  at 
least  one  of  their  nearest  eight  neighboring  points). 
Following  Richtmyer's  suggestion  (I96I),  we  divide  irregular 
net  points  near  the  free  boundary  (front)  into  two  classes. 
A  type  "A"  point  is  an  irregular  net  point  lying  closer  to 
the  boundary  than  a  certain  distance  5,  which  may  be  of  the 
order  of  one  quarter  of  the  grid  distance  A|  (or  Ari);  a  type 
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"B"  point  Is  an  irregular  net  point  not  of  type  A  (see 
Figure  H-)  .      The  procedure  assumes  that  values  of  a 
function  f  =  f{i,T])    are  given  at  regular  net  points,  at 
type  B  points,  and  at  points  on  the  front,  i.e.  every- 
where except  at  type  A  points .   Type  A  points  lie  too 
close  to  the  front,  and  therefore  the  difference  equations 
are  not  used  at  such  points.   Instead,  we  find  the  values 
of  f  at  type  A  points  by  using  quadratic  interpolation. 
For  the  A  point  as  marked  Ap  in  Figure  h,    the 
ordinate  of  intersection  of  the  line  4  -    const, 
through  the  point  Ap  with  the  front  (designated  Y 
in  the  figure)  is  found  by  cubic  interpolation 
along  the  front  knowing  the  coordinates  of  the 
four  consecutive  points  i{-l,  &,    £+1   and  &+2.      A 
value  of  f  at  the  point  Y  is  then  interpolated  by 
cubic  interpolation  along  the  front  using  f  values 
at  the  same  four  consecutive  front  points.   On  the 
line  4  =  const,  through  the  point  Ap,  wa  have  f 
values  at  three  points,  one  regular  and  one  type  B 
point  above  the  point  A^  and  the  Y  point.   We  can 
therefore  obtain  an  f  value  at  the  point  Ap  by 
quadratic  interpolation.   Note  that  we  can  obtain 
another  value  of  f  at  this  point  by  a  similar 
interpolation  procedure  along  the  line  r|  =  const, 
through  the  point  Ap.   As  the  final  value  of  f  at 
the  point  Ap,  we  take  the  arithmetic  mean  of  the 
above  two  f  values.   (If  the  above  quadratic 
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interpolation  is  possible  in  only  one  direction, 

then  we  use  this  value  for  f.) 

2     2 
¥e  are  now  able  to  evaluate  bf/bi,    b   f/S|  ,    5f/Sr|, 

2     2 
B  f/Sri   at  type  A  and  B  points. 

Whenever  net  points  (including  type  A  points) 
on  opposite  sides  of  the  type  B  point  are  available, 
an  ordinary  centered  difference  quotient  is  used. 
Otherwise,  a  point  similar  to  the  point  Y  on  the 
boundary  must  be  used  in  place  of  one  of  the 
opposite  neighbors,  and  then  the  approximation  to 
the  derivatives  must  take  the  unequal  spacing  of 
the  points  into  account,  to  achieve  overall  second- 
order  accuracy. 

The  evaluation  of  S  f/S^Sr|  is  now  carried  out,  at  the 
B  points  only,  with  the  use  of  the  above  calculated  5f/Sr) 
and  Bf/S|  as  follows: 

We  check  whether  the  B  point  has  two  neighbors  in  the 
^  direction  (or  the  t]  direction)  .   If  so,  then  we 
calculate  "1^  (-^)  (o^  ^  ^"sf  ^  ^  ^^  ^  centered  fashionj 
if  not,  then  we  calculate-^  (-^)  (o^-^if:  ("^^^  ^^  ^ 
noncentered  way.   In  other  words,  we  check  through 
the  above  list  of  alternative  ways  of  approximating 
S^f/S^Sri,  and  pick  the  first  method  that  is  possible. 
This  depends  on  the  location  of  the  front.   The  use 
of  an  uncentered  difference  quotient  in  this  case  is 
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permissible  because  the  error  Incurred  In  (6.2)  Is 
of  third  order  since  the  mixed  derivative  terms  are 
multiplied  by  (At)   In  equation  (6.2). 

Now  we  discuss  our  method  of  evaluating  the  gradient 
of  h  at  the  front  point  i    (see  Figure  5) . 

First,  the  direction  of  Vh  at  the  point  I,    Is  determined 
by  drawing  the  Inward  normal  unit  vector  n  at  the 
point,  calculated  by  fitting  a  parabola  to  the  points 
[i-l] ,    £    and  (^+1).   Let  Z  be  the  first  intersection 
of  the  normal  vector  n  and  a  coordinate  line  which 
connects  two  net  points  D..   Also,  let  S  be  the  ('■- 
distance  between  the  point  £   and  the  intersection  Z 
(see  the  inset  of  Figure  5) .   In  general,  one  may 
encounter  the  situation  in  which  the  distance  S  is 
less  than  a  certain  small  distance  5  which  is  of 
the  order  of  one  quarter  of  the  grid  distance.   If 
this  case  occurs,  one  should  choose  as  Z  the'  next 
intersection  of  the  normal  vector  with  another 
coordinate  line.   Our  task  is  to  evaluate  Sh/Sn  at 
the  point  £    (I.e.  Sh/Sn L )  with  second  order  accuracy. 
We  fit  a  parabola  to  the  height,  h  =  hp(n)  where  n  is 
distance  measured  from  C  along  the  normal  vector  n. 
We  then  set  Sh/^nL  =  dhp(0)/dn.   The  parabola  is 
determined  from  three  pieces  of  information,  i.e. 
(dh/dn)  and  h  at  the  point  Z  and  h  =  0  at  the  point  £. 
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The  value  of  "n  at  Z  is  found  by  quadratic  Interpola- 
tion using  the  two  net  points  In  D.  (points  B-,    and  Bp 
In  Figure  5)  on  the  coordinate  line  through  Z  and  a 
third  point  which  may  be  a  regular  net  point  (as  Is 
P  In  the  case  of  Figure  5),  a  type  A  or  B  point,  or 
even  on  the  boundary.   (E.g.,  in  Figure  5  the  third 
point  is  chosen  to  be  P  unless  "ZBp  <  ZB^  in  which 
case  the  third  point  would  be  X.)   The  value  of  dh/dn 
at  Z  is  found  as  follows:   First,  Vh  at  Z  is  obtained 
by  linear  interpolation  using  values  of  Vh  (obtained 
by  the  method  described  in  the  previous  paragraph) 
at  the  two  neighboring  points  B-,  and  B^.   Secondly, 
we  project  Vh  at  Z  onto  the  vector  n,  i.e.,  dh/dn  = 
n-Vh. 

Now,  we  discuss  the  special  treatment  at  the  northern 
boundary.   There,  the  T]-com.ponent  v  of  the  velocity  vanishes 
for  all  time,  but  the  ^-component  u  of  the  velocity  and  the 
height  h  are  computed . 

The  first  and  second  derivatives  of  any  grid  function 
f  with  respect  to  ^  along  the  northern  boundary  are 
evaluated  by  using  schemes  (6.3).   In  order  to 
evaluate  ^r/hr\   and  ^^f/^T)^  so  that  (6.2)  will  be 
accurate  to  second  order  at  a  grid  point  on  the 
northern  boundary,  we  fit  a  parabola  to  the  values 
of  f  at  three  grid  points  nearest  to  the  northern 
boundary  including  one  at  the  boundary.   By 
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differentiating  with  respect  to  t^  the  quadratic 
polynomial  of  f  thus  obtained,  and  evaluating  the 

resulting  derivatives  at  the  point  in  question,  we 

2    2 
obtain  the  evaluations  of  Sf/Sr)  and  ^  f/^r)   along 

the  northern  boundary.   A  difference  scheme  centered 

in  i    is  used  to  compute  S  f/SriS^  =  -rr    (-r-)  at  grid 

points  on  the  boundary,  once  Sf/3r)  has  been  evaluated 

as  above  at  each  grid  point  along  the  boundary. 

Lastly,  a  remark  should  be  made  concerning  the  change 
in  the  set  D.  with  respect  to  time.   The  solution  is 
advanced  in  time  with  (6.2)  at  regular  points  and  at  type  B 
points,  but  not  at  type  A  points  in  the  domain  D  .   The 
movement  of  the  front  points  is  calculated  from  (7-1)  and 
(T-'^)  with  the  steps  described  by  (7-5)  •   At  the  end  of 
step  3  (i.e.,  after  the  first  estimate  of  the  position  of 
the  front  points),  we  check  for  any  change  in  the  set  D. . 
This  is  done  as  follows: 

First  compute  the  locations  of  the  intersection  of 
the  front  with  the  ^,t]  coordinate  lines,  and 
secondly  inspect  the  domain  index  of  the  grid  points 
in  the  new  domain  D.  at  time  t  +  Ax.   For  this  purpose, 

we  use  a  zero  index  to  identify  grid  points  in  the 

t 

warm  air  domain  D.    and  a  nonzero  index  to  denote  the 

set  in  the  domain  D. .   At  the  end  of  step  3,  after 
checking  through  the  indices  of  points  in  the  new 
domain  D  at  time  t  +  Ax,  if  one  or  more  zero  indexed 
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points  occur,  then  they  must  be  relabeled  as  nonzero 
Indexed  points  and  classified  as  type  A  points.   In 
the  case  that  the  domain  D,  is  locally  expanding,  a 
type  B  point  at  time  x   m.ay  be  eligible  as  a  regular 
point  at  t+At,  similarly  a  type  A  point  may  change 
to  a  type  B  point.   In  the  case  that  the  domain  D.  is 
locally  contracting,  a  regular  point  at  time  x  may 
degenerate  to  a  type  B  point  at  t+At  and  similarly  a 
type  B  point  may  change  to  a  type  A  point.   This 
procedure  presupposes  that  At/As  is  sufficiently  small 
so  that  during  the  movement  of  the  front  from  time  t 
to  t+At,  we  gain  or  lose  only  type  A  points  when  any 
changes  occur  in  the  set  D. . 

8.   Data  and  results  of  the  calculations. 

The  following  numerical  values  for  the  parameters  in 
the  problem  were  taken. 

As  =  250,000  feet  (=  76.2  ki-n) 

At  =   600  sec,  A  =  .0024  sec/feet 

X  =   20  As 

Y  =   20  As 

C^   =   2  As 

Cg  =   9.5  As 

g  =  32.1521  feet/sec^  {=   9.8m/sec^) 

-4    -1 
f  =  10   sec 

u  =   10  feet/sec 

£—  u   =50  feet/sec 

g(l  -  ^)      =   0.6  feet/sec^ 
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The  above  magnitude  of  the  density  discontinuity 
corresponds  to  a  temperature  discontinuity  of  the  order 
of  5  degrees  Centigrade.   The  slope  Sh/Sy  of  the 
stationary  discontinuity  surface  as  given  by  the  formula 
(2.5b)  Is  found  to  have  the  value  I/I50  when  the 
parameters  are  chosen  as  above,  and  that  Is  In  the  range 
of  observed  values  for  the  slope  of  frontal  surfaces  In 
middle  latitudes.   The  time  step  At  -    6OO  sec  was  chosen 
on  the  basis  of  the  convergence  criterion  described  In 
Section  6. 

We  shall  first  describe  the  results  of  calculation 
for  Case  A.   Figure  6a  shows  the  Initial  height  contour 
pattern  of  the  cold  air.   The  contour  lines  h  =  constant 
are  drawn  at  5,000  foot  Intervals.   The  front,  given  by 
a  sinusoidal  curve,  separates  the  domain  of  the  cold  air 
from  that  of  the  warm  air.   The  solid  circles  on  the 
front  Indicate  the  Initial  locations  of  the  front  points 
and  are  uniformly  spaced  at  Intervals  Ax  =  0.8  As.   The 
marks  along  the  boundaries  of  the  entire  domain  R  and  the 
Inset  of  Figure  6a  show.  In  correct  scale,  the  size  of 
the  grid  lattice.   For  convenience  In  programming,  two 
extra  columns  of  grid  points  were  added  outside  of  the 
western  boundary  and  one  extra  column  of  grid  points 
outside  of  the  eastern  boundary  to  handle  the  periodicity 
conditions  for  the  dependent  variables  at  those  boundaries . 

Figure  6b  shows  the  contour  pattern  of  the  cold  air 
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height  at  t  =  8  hours,  when  observed  from  the  "moving 
coordinate"  system,  i.e.  as  seen  from  a  coording,te  system 
translated  8  x  360O  x  10  feet  to  the  east  of  the  initial 
position.   (All  the  diagrams  presented  in  this  section 
refer  to  the  moving  coordinates  discussed  in  Section  5.) 
The  location   of  the  initial  curve  of  the  front  corres- 
ponding to  a  uniform  translation  of  10  feet/sec  is 
indicated  by  the  dashed  line  connected  with  open  circles. 
Evidently  the  entire  front  progressed  eastward  relative 
to  the  moving  coordinates .  Also  the  cold  front  moved 
eastward  faster  than  the  warm  front.   The  development  of 
this  characteristic  agymmetry  suggests  the  occlusion 
process  of  frontal  cyclones  (BJerknes,  J.  and  H.  Solberg, 
1922),  which  agrees  with  the  results  of  qualitative 
analyses  given  by  '    Whitham  (1955)  and  Stoker  (1955,7). 

In  order  to  show  the  movement  of  the  front  in  detail, 
the  trajectories  of  individual  points  on  the  front  during 
the  period  of  8  hours  are  shown  in  Figure  7.   Note  that 
points  on  the  cold  front  moved  southeastward  and  those  on 
the  warm  front  moved  northwestward  on  the  average,  whereas 
both  fronts  themselves  propagated  eastwards.   The  movement 
of  points  on  the  front  clearly  indicates  the  production  of 
a  cyclonic  circulation  around  the  circulation  center  at 
which  the  cold  front  joins  the  warm  front.   The  magnitude 
of  the  velocity  components  (in  the  coordinate  system  moving 
eastward  with  velocity  10  feet/sec)  of  the  front  points  at 
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t  =  8  hours  are  shown  by  the  numerals  in  units  of 
feet/sec.   The  numerators  are  the  x-component  of 
velocity  and  the  denominators  the  y-component. 

The  trajectories  of  the  points  of  the  front  near 
the  circulation  center  show  the  converging  motion  of 
the  cold  air.   Consequently,  the  spacing  between  these 
points  becomes  small  compared  with  the  grid  spacing. 
In  part  owing  to  this  uneven  distribution  of  these 
points  on  the  front  and  in  part  owing  to  a  sharp 
transition  in  the  orientation  of  the  front  near  the 
circulation  center,  the  cubic  interpolation  procedure 
to  approximate  the  curve  of  the  front  becomes  inaccurate, 
and  so  do  the  trajectory  calculations.   A  recipe  for 
avoiding  this  congestion  of  points  on  the  front  is  to 
redefine  their  location  from  time  to  time  so  that  two 
consecutive  points  are  neither  too  close  together  nor  too 
far  apart.   In  the  two  problems  solved  numerically  here 
no  redistribution  of  points  on  the  front  was  made,  but 
this  technique  worked  out  successfully  in  other  cases  of 
the  same  kind.   It  was  found,  however,  that  ultimately 
difficulties  arise  not  only  from  the  uneven  spacing  of 
front  points,  but  also  from  a  loss  in  accuracy  in  approxi- 
mating the  curve  of  the  front  by  cubic  interpolation  in 
the  region  of  high  curvature.   A  remedy  for  this  latter 
inaccuracy  would  be  the  introduction  of  a  finer  spacing 
of  points  along  the  front,  at  the  expense  of  reducing  the 
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time  increment  of  computation;  such  a  procedure  is  planned 
for  the  future . 

The  results  of  calculation  of  Case  B  are  now 
described.   Figure  8a  shows  the  Initial  height  contour 
pattern  of  the  cold  air,  this  is  similar  to  Figure  6a. 
Since  it  is  desired  to  have  the  initial  wind  field 
geostrophic  (see  (j.Jh))  and  the  northern  boundary 
condition  is  that  v  =  0,  the  height  along  the  northern 
boundary  is  constant  initially.   The  slope  of  the  height 
of  cold  air  Sh/Si]  is  constant  with  respect  to  t],  but  it 
is  variable  with  respect  to  ^. 

Figure  8b  shows  the  height  contour  pattern  of  the 
cold  air  at  t  =  11  hours.   Just  as  in  Case  A  (see  Figure 
6b),  it  is  again  observed  that  the  entire  front  system 
progressed  eastward  relative  to  the  moving  coordinate 
system.   The  configuration  of  the  front  clearly  indicates 
a  development  tending  toward  occlusion. 

Figure  9  illustrates  the  motion  of  the  front  by 
showing  the  trajectories  of  individual  front  points.   As 
in  Case  A,  the  circulation  center  progresses  eastward 
relative  to  the  moving  coordinate  system  and  overtakes 
the  warm  front  leaving  warm  air  behind  the  cold  front. 
As  discussed  previously,  the  trajectories  of  the  points 
on  the  front  near  the  circulation  center  indicate  not 
only  the  production  of  a  cyclonic  circulation,  but  also 
the  converging  motion  of  cold  air.   Converging  air  must 

-  33  - 


be  lifted  upwards,  and  the  subsequent  thermodynamic 
processes  (such  as  condensation  of  water  vapor  by 
adlabatlc  cooling)  produce  moisture.   Hence  the  presence 
of  a  relatively  strong  mass  convergence  behind  the  cold 
front  Implies  that  severe  storms  are  likely  to  be 
associated  with  cold  fronts. 

As  one  check  of  the  numerical  calculations,  the 
magnitudes  of  the  total  mass  M  defined  In  (4.1),  the 
total  kinetic  energy  K  and  potential  energy  P  of  the 
cold  air  defined  In  connection  with  (4.2),  and  the  rate 
of  total  work  W,  done  upon  the  cold  air  layer  by  the 
action  of  the  Corlolls  force  In  the  warm  air  layer, 
defined  by  (4.3)  were  all  computed  at  every  time  step. 
For  both  Cases  A  and  B,  the  total  mass  M  was  conserved 
very  well;  the  variations  amounted  to  less  than  0.01/6. 
For  both  Cases  A  and  B,  the  values  of  W  stay  negligibly 
small  during  the  period  of  calculation.   This  circumstance 
Is  due  mainly  to  our  choice  In  the  Initial  conditions. 
Thus,  the  sum  of  the  potential  and  kinetic  energy  remained 
approximately  constant  with  the  variation  0.02/0  In  case  A 
and  0.2/o  In  Case  B.   (The  value  of  W  In  Case  B  was  roughly 
one  order  of  magnitude  larger  than  that  In  Case  A.) 

In  Figure  10  the  values  of  K  In  Cases  A  and  B  were 
plotted  In  units  of  p (As/At)   x  10"   as  a  function  of  time. 
Here  the  total  kinetic  energy  was  computed  by  using  the 
velocity  of  the  wind  measured  In  the  moving  coordinate  system. 

-  34  - 


In  Case  A,  the  relative  wind  field  was  zero  Initially, 
so  that  the  value  of  K  at  t  =  0  Is  zero.   The  values  of 
K  Increased  steadily  up  to  about  t  =  4  hours  and  then 
leveled  off  after  that.   The  source  of  kinetic  energy 
is  of  course  the  release  of  the  potential  energy  of  the 
cold  air.   On  the  other  hand,  in  Case  B  the  calculation 
started  from  a  state  of  geostrophic  balance.   The 
exchange  of  energy  between  the  kinetic  and  potential 
forms  appears  to  be  nearly  periodic  with  the  period 
of  about  10  hours . 

9 .   Conclusions . 

This  article  represents  a  first  step  in  the 
development  of  numerical  methods  for  the  study  of 
frontal  cyclones  in  the  atmosphere.   In  this  attempt, 
to  simplify  the  problem,  the  motion  of  the  lower  cold 
air  layer  has  been  treated  by  assuming  that  the  dynamics 
of  the  perturbations  in  the  upper  warm  layer  could  be 
neglected  (see.  Section  2). 

The  two  cases  studied  in  this  work  demonstrate 
that  the  action  of  the  purely  mechanical  forces,  gravity 
and  Coriolis,  are  primarily  responsible  for  the  gross 
features  of  the  occlusion  process. 

These  studies  are  being  continued  by  adding  the 
warm  air  layer  to  the  dynamical  system  and  considering 
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the  coupling  of  the  two  layers.   This  Is  done  by  Integrat- 
ing numerically  the  differential  equations  for  such  a 
system  presented  In  Section  2.   This  model  may  throw  some 
light  on  the  nature  of  the  Interaction  between  the  low 
tropospherlc  frontal  cyclones  at  the  polar  front  and  the 
upper  long  waves  or  planetary  waves  (as  defined  by  Rossby 
and  collaborators  (1939))'   The  nature  of  the  interaction 
has  been  investigated  synoptlcally  by  Palmen  (1951), 
Berggren,  Bolln  and  Rossby  (1951),  Palmen  and  Newton  (I95I) 
and  others.   This  has  been  summed  up  by  J.  Bjerknes  (I95I) 
as  follows:   "we  may  classify  the  formation  of  extra- 
trpoical  cyclones  as  being  due  either  to  unstable  frontal 
wave  action  or  to  unstable  growth  of  an  upper  wave  trough. 
A  subsequent  combination  of  both  processes  is  quite 
frequentj  and  all  the  strongest  cyclones  on  record  seem 
to  have  that  double  origin." 


-  36 


Acknowledgement 

We  are  very  grateful  to  our  colleagues  at  the  Courant 
Institute  of  Mathematical  Sciences,  New  York  University, 
for  their  generous  advice  and  assistance  in  completing  this 
study.   In  particular,  we  received  many  useful  comments  in 
regard  to  the  formulation  of  the  problem,  finite-difference 
schemes  and  machine  calculations  from  George  K.  Morikawa, 
Robert  D.  Richtmyer,  Peter  D.  Lax  and  Max  Goldstein.   Most 
of  the  numerical  computations  were  performed  on  the 
IBM  7090  com.puter  at  the  AEC  Computing  and  Applied 
Mathematics  Center  of  the  Institute.   Eva  Swenson  assisted 
us  in  programjning  a  part  of  the  code. 

This  manuscipt  was  completed  after  one  of  the  authors 
(Kasahara)  moved  to  the  National  Center  for  Atmospheric 
Research,  Boulder,  Colorado.   Further  num.erlcal  calcula- 
tions were  carried  out  on  the  CDC  360O  at  the  Computing 
Center  there.   We  record  here  our  appreciation  to 
Dr.  Glenn  E.  Lewis,  Director  of  the  Computing  Center  (who 
had  developed  an  earlier  version  of  the  difference  scheme 
for  the  cold  front  model,  mentioned  in  the  Introduction, 
when  he  was  at  the  Courant  Institute).   In  programming 
we  were  aided  by  Charmian  Jefferies.   The  authors  take 
pleasure  in  acknowledging  the  cooperation  given  us  by 
the  National  Center  for  Atmospheric  Research. 


37 


Appendix 

The  form  of  S  and  K  In  (6.3)- 

Let  us  Introduce  abbreviations  such  as 


P^     =  P(|,i1,t)  ;   P'^"''-'"  -  ?{i}y],x+A'v)      etc 


The  equation  (6.2)  is  then  written  as 

p'^+1   -   P""  -I  (A^^^  +  A^IJ^P-^  AT 
(A.l)  -i    (b'^+I  +  b'^)^P^  At 


+  -|  Q(P^''"-^  +  P'^)  At  +  KAt 


where 

J;  =    a      At  ,    _af     aA  a  5^    ,  sb  a      ^  a  ^ 

li  =  ^   At  ,    5^   ,  ^A  a  ,  „   a^  ,  as  a    „  a  ^ 

^  -  ^  -  ~  ^^^T^  "^^^  +  ^^ +^^  -  Q^^ 

By  transposing  quantities  at  time  level  t+At  to  the  left 
hand  side  of  (A.l)  and  quantities  at  t  to  the  right  hand 
side,  we  obtain 

p-^+l  +1  (A^+I^P-^  +  B^+l^P^  -  Q  p-^+l)  AT 
(A. 2) 

=     {  I  -  -i  (a'^  ^  +  b'^  5/  -  Q)  J  P'^  AT  +  KAT  . 

Since  the  elements  of  A  and  B  are  linear  functions  of  the 
components  of  P,  we  may  define  the  matrix  M  such  that 
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(A. 3)     M  P^+1  =  AE^A^+i^p-^  +  b^+Vp''  -  Q  P''"'^] 


with  the  aid  of  (A. 3),  (A. 2)  may  be  expressed  as 


P^+^  =   S  P^  AT  +  K*  AT 


where 


S   =   (I  +  M)"^  [I  -  ^  (a"^-^  +  ^^y-    -   Q)At] 


K*  =   (I  +  M)"^  K 


and 


M  =  ^ 
2 


/^u^ 


47v'^+P 


y 
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Legend 

Pig.  la    A  stationary  front. 

Fig.  lb    Vertical  cross  section  of  the  two  layers. 

Fig.  2     Flows  in  cold  and  warm  air  layers. 

Fig.  3     The  rectangular  domain  of  numerical 

integration. 
Fig.  4     Classification  of  the  type  of  grid  point 

near  the  frontal  boundary. 
Fig.  5     Evaluation  of  Vh  at  the  frontal  boundary. 
Fig.  6a    Height  contour  pattern  of  cold  air  for 

Case  A  at  t  =  0 .   The  contour  lines  are 

drawn  at  5,000  feet  intervals. 
Pig.  6b    Height  contour  pattern  of  cold  air  for 

Case  A  at  t  =  8  hours.   The  contour  lines 

are  drawn  at  5,000  feet  intervals. 
Pig.  7     Trajectories  of  front  points  (portion) 

during  the  period  of  8  hours  for  Case  A. 
Pig.  8a    Height  contour  pattern  of  cold  air  for 

Case  B  at  t  =  0,  similar  to  Pig.  6a. 
Fig.  8b    Height  contour  pattern  of  cold  air  for 

Case  B  at  t  =  11  hours,  similar  to  Fig.  6b. 
Fig.  9      Trajectories  of  front  points  during  the  period 

of  11  hours  for  Case  B.   Numerals  beside  each 

front  point  are  the  identification  numbers. 

Fig.  10    Total  kinetic  energy  as  a  function  of  time. 
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Pig.  lb.  Vertical  cross-section  of  the  two  layers. 
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Fig.  2.   Flows  In  cold  and  warm  air  layers. 


<0 

< 

Ul 


^  k 


Q 


O 


■Q 


o 


< 
o 

O 
O 


0£ 

< 


01 
< 


FRONT 

D  TYPE  A 

A  TYPE  B 

O  REGULAR  NET  POINTS 

•  FRONT  POINTS 


IRREGULAR  POINTS 


Pig.  4.   Classification  of  the  type  of  grid  point  near 
the  frontal  boundary. 
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Fig.  6a.  Height  contour  pattern  of  cold  air  for  Case  A  at  t  =  0 
The  contour  lines  are  drawn  at  5,000  feet  Intervals. 
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'ig .  6b.  Height  contour  pattern  of' cold  air  for  Case  A  at  t  =  8  hours 
The  contour  lines  are  drawn  at  5,000  feet  Intervals. 
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Fig.  8a.  Height  contour  pattern  of  cold  air  for  Case  B  at  t  =  0,  similar 
to  Pig.  6a. 
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